#include <stdio.h>
#include <math.h>
#include <sys/time.h>
#include <omp.h>
#include "hw2.h"
#include <stdlib.h>
#include <string.h>

#define pi 3.1415926535877932384626
#define G 0.001

//Wall clock time in seconds
double
wctime ()
{
  struct timeval tv;
  gettimeofday (&tv, NULL);
  return (tv.tv_sec + 1E-6 * tv.tv_usec);
}

/** initialize mass, velocity and position */
void
initial (double m[], double v[][2], double p[][2], int n)
{
  int i;
  for (i = 0; i < n; i++)
    {
      m[i] = 1;
      double theta = i * 2 * pi / n;
      // distributed around unit circle
      p[i][0] = cos (theta);
      p[i][1] = sin (theta);
      // velocity point to origin, or normal direction
      v[i][0] = -p[i][0];
      v[i][1] = -p[i][1];
    }
}

int
main (int argc, char **argv)
{
  double m[N];			//Mass
  double v[N][2];		//Velocity
  double p[N][2];		//Position
  double ftot[N][2];		//Total force
  double a[N][2];

  int n;			//Number of bodies

  int i, j, k;
  double dt = 0.0001;		//Time interval
  double mv0 = 0, mv1 = 0, mvr = 0;

  double iter = 4;

  double begin, end;
  double timecost;

  if (argc > 2)
    {
      n = atoi (argv[1]);
      if (n > N)
	{
	  printf ("Error: the input n is too large!\n");
	  return 1;
	}
    }
  else
    {
      printf ("Error: must indicate body number and output file name!\n");
      return 1;
    }

  int num_threads = omp_get_max_threads ();

  double (*thread_accumu)[N * 2];
  thread_accumu = malloc (sizeof (double) * N * num_threads * 2);

////////////////////////////////////////////////////////////////////////


  initial (m, v, p, n);

  //Compute the net momentum and angular momentum
  for (i = 0; i < n; i++)
    {
      mv0 += m[i] * v[i][0];
      mv1 += m[i] * v[i][1];
      mvr += m[i] * (p[i][0] * v[i][1] - p[i][1] * v[i][0]);
    }

  printf ("Pre :\tmomentum: %f  %f\n\tanglular momentum: %f\n", mv0, mv1,
	  mvr);

  begin = wctime ();

  for (k = 0; k < iter; k++)
    {

      memset (ftot, 0, sizeof (double) * N * 2);
      memset (thread_accumu, 0, sizeof (double) * N * num_threads * 2);

#pragma omp parallel
      {
	//compute f_ij for all elements
// #pragma omp for schedule(static) private(i, j)
//      for (i = 0; i < n; i++)
//        {
//          for (j = 0; j < num_threads; j++)
//            {
//              thread_accumu[j][2 * i] = thread_accumu[j][2 * i + 1] = 0;
//            }
//          ftot[i][0] = ftot[i][1] = 0;
//        }

#pragma omp for private(i, j)  schedule(guided)
	for (i = 0; i < n; i++)
	  {
	    int thread_no = omp_get_thread_num ();

	    for (j = i + 1; j < n; j++)
	      {
		double r_ij0, r_ij1, r_ij, rcube, f_ij0, f_ij1, f_ij;
		r_ij0 = p[j][0] - p[i][0];
		r_ij1 = p[j][1] - p[i][1];
		r_ij = sqrt (r_ij0 * r_ij0 + r_ij1 * r_ij1);
		rcube = r_ij * r_ij * r_ij;
		f_ij0 = G * m[i] * m[j] * r_ij0 / rcube;
		f_ij1 = G * m[i] * m[j] * r_ij1 / rcube;

		ftot[i][0] += f_ij0;
		ftot[i][1] += f_ij1;

		thread_accumu[thread_no][j * 2] -= f_ij0;
		thread_accumu[thread_no][j * 2 + 1] -= f_ij1;
	      }
	  }

#pragma omp for private(i, j) schedule(static)
	for (i = 0; i < n; i++)
	  {
	    for (j = 0; j < num_threads; j++)
	      {
		ftot[i][0] += thread_accumu[j][i * 2];
		ftot[i][1] += thread_accumu[j][i * 2 + 1];
	      }
	  }

	//update velocity and position
#pragma omp for schedule(static) private(i)
	for (i = 0; i < n; i++)
	  {
	    a[i][0] = ftot[i][0] / m[i];
	    a[i][1] = ftot[i][1] / m[i];

	    p[i][0] += v[i][0] * dt;
	    p[i][1] += v[i][1] * dt;

	    v[i][0] += a[i][0] * dt;
	    v[i][1] += a[i][1] * dt;
	  }

      }				//End Parallel region

    }				//End iteration

  end = wctime ();

  free (thread_accumu);

  //Compute the net momentum and angular momentum
  mv0 = mv1 = mvr = 0;
  for (i = 0; i < n; i++)
    {
      mv0 += m[i] * v[i][0];
      mv1 += m[i] * v[i][1];
      mvr += m[i] * (p[i][0] * v[i][1] - p[i][1] * v[i][0]);
    }

  printf ("Post:\tmomentum: %f  %f\n\tanglular momentum: %f\n", mv0, mv1,
	  mvr);

  timecost = end - begin;
  double t_iter = timecost / iter;	//Time per iteration

  FILE *fp;
  fp = fopen (argv[2], "a");
#ifdef _OPENMP
  fprintf (fp, "%d ", num_threads);
#endif
  fprintf (fp, "%d ", n);
  //fprintf (fp, "Time: %f\n", timecost);
  fprintf (fp, "%f\n", t_iter);
#ifdef _OPENMP
  printf ("Num threads: %d\n", num_threads);
#endif
  printf ("Body num: %d\n", n);
  printf ("Time: %f\n", timecost);
  printf ("Time per iteration: %f\n\n", t_iter);
  fclose (fp);

  return 1;

}
